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INVESTIGATION  OF  AN  ALTERNATIVE  FINITE  ELEMENT  PROCEDURE: 
A  ONE-STEP,  STEADY-STATE  ANALYSIS 


Leonard  R.  Herrmann  and  Joe  Mello 
Department  of  Civil  and  Environmental  Engineering 
University  of  California 
Davis,  CA.  95616 


Introduction 

The  purpose  of  this  study  was  to  investigate  the  feasibility  of  developing  a 
one-step,  steady-state  finite  element  analysis  procedure  that  is  applicable  to 
the  "plow  problem".  The  plow  problem  being  defined  as  the  determination 
of  the  state  of  stress  and  deformation  induced  in  a  soil  deposit  by  the  passage 
of  a  plow  type  device  through  the  soil  at  a  constant  velocity  and  the 
calculation  of  the  force  required  to  drive  the  plow.  The  analysis  of  this 
problem  can  be  carried  out  as  a  transient  analysis  using  available  commercial 
finite  element  codes.  Such  an  analysis  would,  in  a  step  wise  fashion,  analyze 
the  configuration  as  the  plow  entered  the  soil  and  proceeded  to  reach  a  steady- 
state  plowing  action.  The  advantage  of  such  an  analysis  is  that  it  can  be 
performed  using  available  software.  In  addition,  if  the  plow  insertion  is 
described  in  a  realistic  fashion,  the  stresses  and  deformations  developed 
during  the  insertion  process  would  also  be  determined.  The  only  possible 
drawback  to  such  an  analysis  is  possible  large  computational  costs  associated 
with  a  transient,  nonlinear,  inelastic,  three-dimensional  finite  element 
analysis  of  a  relatively  complicated  configuration. 

Assuming  that  the  steady-state  plowing  process  does  not  result  in  any 
periodic  localization  (e.g.,  cracks  propagating  into  the  soil  mass  at  regular 
intervals  from  the  plow  path),  an  alternative  would  be  to  perform  a  one-step, 
steady-state  analysis.  An  observer  situated  on  a  plow  passing  through  a 
homogeneous  soil  mass  at  a  constant  velocity  would  observe  steady-state 
conditions.  The  goal  of  a  steady-state  analysis  is  to  capture  this  steady-state 
behavior  in  a  one  step  analysis.  The  advantage  of  such  an  analysis  is  that 
multiple  solution  steps  are  not  required  and,  thus,  there  is  a  potential  for 
considerable  computational  cost  saving.  The  disadvantage  is  that  when  the 
problem  is  nonlinear  due  either  to  material  or  geometric  nonlinearities, 
commercial  software  does  not  exist  to  perform  the  analysis.  Approximate 
analyses  [1,2]  for  specialized  geometric  configurations  have  been  reported  for  a 
related  class  of  Geotechnical  Engineering  problems  (insertion  of  a  sampling 
tube  into  a  soil  mass),  however,  in  order  to  determine  the  deformation 
pattern  the  method  treats  the  soil  as  a  liquid;  the  degree  of  approximation 
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introduced  by  this  assumption  is  hard  to  access.  In,  addition,  no  general 
procedures  are  available  to  extended  this  type  of  analysis  to  the  plow  problem. 

Weighting  the  advantages  and  disadvantages  of  the  two  competing  analysis 
procedures  it  seemed  clear  that  the  available  and  proven  multi-step,  transient 
analysis  procedure  is  the  preferable  way  to  proceed  for  a  project  requiring  a 
solution  within  a  given  time  frame.  However,  as  a  possible  backup  in  case 
the  transient  analysis  should  prove  to  be  excessively  expensive,  it  seemed 
advisable  to  perform  a  preliminary  investigation  of  the  feasibility  of  the 
steady-state  procedure.  Such  a  preliminary  investigation  was  the  goal  of  this 
portion  of  the  project. 


Scope  of  Study 

In  order  to  model  the  plowing  process  as  steady  state  the  following 
assumptions  (these  may  be  viewed  as  restrictions  or  approximations)  are 
required; 

1)  It  is  assumed  that  the  soil  deposit  is  infinite  in  extent,  level  and 
homogeneous  in  the  direction  of  the  plowing  action. 

2)  It  is  assumed  that  the  plow  moves  at  a  constant  depth,  with  constant 
velocity  Vq  and  in  a  straight  path. 

3)  It  is  assumed  that  no  significant  "periodic  localizations"  occur  in  the 
soil  mass  during  the  plowing  process.  (If  significant  "periodic 
localizations",  should  occur,  then  an  observer  on  the  plow  would  not 
observe  steady-state  conditions.)  An  example  of  a  periodic  localization 
would  be  the  radiation,  from  the  plow  path,  of  cracks  into  the 
surrounding  soil  at  regular  intervals  along  the  path. 

The  objective  of  this  study  was  to  demonstrate  the  feasibility  of  the  one-step, 
steady-state  analysis  procedure  for  the  plow  problem.  In  order  to  carry  out  the 
investigation  within  the  time  and  financial  constraints  of  the  project, 
additional  restrictions  were  introduced  in  order  to  simplify  the  analysis. 
However,  it  should  be  emphasized  that  these  restrictions  are  not  required  for 
a  steady-state  analysis  to  apply.  It  is  hoped  that  even  within  the  confines  of 
these  simplifying  assumptions  that  all  the  basic  features  of  a  one-step,  steady- 
state  analysis  can  be  demonstrated  and  that  success  for  this  restricted  class  of 
problems  will  indicate  the  likelihood  of  success  if  the  method  were  to  be 
applied  to  the  actual  plow  problem. 

The  additional  simplifying  assumptions  are: 
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1)  The  geometry  will  be  modeled  as  2-dimensional  (plane  stress  or 
plane  strain)  rather  than  3-dimensional  as  is  the  case  for  the  actual 
plow  problem. 

2)  It  will  be  assumed  that  the  shape  of  the  plow  is  such  that  the  soil 
remains  in  contact  with  the  plow  along  its  entire  surface  and  that  no 
crack  in  the  soil  propagates  and  opens  ahead  of  the  plow. 

3)  It  is  assumed  that  the  velocity  Vq  is  sufficiently  large  and  the 
permeability  of  the  soil  is  sufficiently  small,  that  the  soil  mass  can  be 
assumed  to  be  "undrained”. 

4)  The  soil  will  be  modeled  as  a  linear  viscoelastic  material 

The  two  neglected  phenomenon  of  pore  water  flow  due  to  the  development 
of  excess  pore  water  pressure  and  the  elastic-plastic  behavior  of  real  soil,  both 
lead  to  a  history  dependency  of  the  solution.  The  assumed  linear 
viscoelasticity  of  the  soil  behavior  also  leads  to  a  history  dependency  of  the 
solution.  It  is  not  intended  to  suggest  that  these  two  history  dependencies  are 
in  any  way  equivalent  (although  one  might  try  to  calibrate  the  viscoelasticity 
model  so  as  to  crudely  capture  the  effects  of  redistribution  of  pore  water 
pressure  and  soil  plasticity).  Instead  it  is  hoped  that  the  demonstration  of  the 
ability  of  the  steady-state  analysis  to  capture  the  history  dependency  of  the 
viscoelastic  properties  will  demonstrate  the  feasibility  of  capturing  the  history 
dependency  of  pore  pressure  redistribution  and  soil  plasticity  by  means  of  a 
steady-state  analysis. 

While  it  is  anticipated  that  inertia  effects  will  only  be  marginally  important 
for  the  plow  speeds  of  interest  in  this  study,  they  are  included  in  order  to 
demonstrate  the  ease  with  which  they  may  be  modeled. 


Definition  of  Problem 

The  purpose  of  this  steady-state  analysis  is  to  study  the  disturbance  in  a  soil 
mass  produced  by  a  plow  and  to  determine  the  force  required  to  drive  the 
plow.  A  simple  two-dimensional  configuration  of  such  a  process  is  shown  in 
Figure  1.  The  design  parameters  are  the  shape  of  the  plow  g,  the  coefficient  of 
friction,  f,  between  the  plow  and  the  soil,  and  the  velocity  of  the  plow  vq- 
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Figure  1.  Two-dimensional  plowing  of  a  soil  mass 


Steady-State  Behavior 

The  initial  insertion  of  the  plow  into  the  soil  is  a  transient  process  which  is 
not  addressed  in  this  study,  i.e.,  it  is  assumed  that  the  plowing  process  has 
proceeded  for  a  long  enough  period  of  time  that  steady-state  conditions  have 
been  reached.  The  analysis  looks  at  the  process  at  a  particular  instant  in  time, 
t=T,  when  the  plow  has  reached  the  point,  relative  to  the  fixed  x-y  coordinate 
system,  shown  in  Figure  2  (an  alternative  point  of  view  is  to  consider  the 


Figure  2.  Plowing  process  at  Time  T 


coordinate  system  as  fixed  to  the  moving  plow).  Point  pc  denotes  the  point  in 
the  soil  that  is  at  the  tip  of  the  plow  at  the  instant  of  time  of  interest  T.  Its 
location  in  space  in  the  undeformed  soil  has  been  selected  to  be  the  origin  of 
the  x-y  coordinates.  The  displacement  the  soil  point  in  question  experiences 
in  going  from  its  undeformed  location  to  its  deformed  location  is  Uc,  see 
Figure  3. 


Figure  3.  Deformed  and  undeformed  grids  at  Time  T 


In  Figure  3  is  shown  a  simple  rectangular  grid  attached  to  the  soil.  Two  views 
are  shown,  one  before  the  soil  is  disturbed  by  the  plow  (undeformed  grid)  and 
the  second  of  the  deformed  soil  at  the  instant  of  time  of  interest,  T.  Consider 
the  generic  node  point  P  fixed  to  the  soil.  With  time,  as  the  plow  moves 

through  the  soil,  the  locations  of  points  P9,  P8,  P7, . (Pi  being  the  last  node 

point  in  the  sequence  at  the  right  extreme  of  the  grid)  will  successively  bear 
the  same  relationship  to  the  location  of  the  plow  as  P  does  currently.  Thus,  if 
one  wishes  to  know  the  past  history  experienced  by  the  soil  at  P  one  merely 
needs  to  look  at  the  present  states  of  points  Pn.  A  similar  statement  may  be 
made  for  the  sequence  Q,  Q9,  Qs,  etc. 

Specifically  the  state  of  the  soil  at  point  P  (located  at  x,y)  at  a  past  time  T-At  is 
identical  to  the  current  state  of  the  soil  at  a  point  located  at  x+voAt,y.  Thus, 
for  example  the  strain  {£}  of  the  soil  at  location  x,y  and  at  time  T-At  is  equal  to 
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the  strain  at  point  x+VoAt,y  at  time  T  (where  {8}  contains  the  strain 
components,  Eij,  written  in  vector  form)  ,  i.e.. 


{e(x,y,T-At)}  =  {e(x+voAt,y,T)} 


(1) 


Thus,  partial  derivatives  with  respect  to  time  are  equal  to  the  partial 
derivatives  with  respect  to  x,  e.g.. 


9u  3u 

aF  =  ■'"o 


(2) 


Governing  Equations 

The  study  will  be  limited  to  small  rotation  and  small  deformation  conditions. 
The  equations  of  motion  are: 


a^xy  _ 

ax  dy  “  P  at^ 

(3) 

axxy  axyy  a^uy 

ax  ay  “  P  at2 

Using  the  results  of  eq  (2)  permits  the  accelerations  to  be  replaced  by  space 
derivatives,  yielding: 


aXxy  aXyy  ^  a^Uy 

-gf  +  -^-pv„2-^=0 

The  deformation  is  described  by  the  strain  terms: 
aux 


£xx  — 


ax 


aUy 

^yy=  ay 


aux  au^ 


Xyv  — 


^y  ay  dx 


(4) 


(5) 
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For  an  inelastic  material,  the  state  of  stress  at  a  particular  point  x,y,  and  at  the 
current  time  T,  is  a  function  of  the  history  of  deformation  at  the  point  in 
question,  i.e., 

{a(x,y,T)}  ={(y({e(x,y,t=-oo  ->  T)})}  (6) 

Where  {a}  is  the  stress  written  in  vector  form.  The  actual  form  of  this 
dependency  depends  upon  the  form  of  the  material  model  used  for  the  soil. 
In  this  study  a  linear  viscoelastic  law  is  used  (e.g.,  see  [3]). 

Using  eq  (1)  the  time  dependence  in  eq  (6)  can  be  expressed  as  a  space 
dependence 

{a(x,y,T)}  ={a({e(x'=x  ->  oo,  y,T)})}  (7) 

The  actual  range  of  x’  must  only  be  carried  far  enough  that  undisturbed  (by 
the  action  of  the  plow)  soil  is  reached. 


Cutting  Condition 

The  soil  is  sliced  to  form  two  new  surfaces  at  the  tip  of  the  plow.  Thus,  in  the 
soil  right  at  the  tip  of  the  plow  some  critical  state  for  new  surface  formation 
must  exist,  this  will  be  referred  to  as  the  "cutting  criterion"  and  must  be 
specified  in  the  finite  element  analysis  of  the  plowing  process.  There  is  a 
problem,  however,  in  that  very  little  appears  to  be  known  about  this  criterion. 
Considerable  work  has  been  done  concerning  the  nature  of  the  failure  state  at 
the  tip  of  an  advancing  free  surface  crack  in  materials  such  as  metals, 
however,  very  little  appears  to  be  known  for  the  corresponding  state  in  soil 
where  the  advancing  crack  tip  is  being  forced  by  a  cutting  tool.  In  the  absence 
of  any  concrete  information,  two  possible  criteria  will  be  considered.  The  first 
is  that  the  normal  effective  stress  across  the  tip  of  the  advancing  cut  is  zero. 
As  an  alternative,  the  second  criterion  is  that  the  soil  must  be  in  a  state  of 
failure,  as  defined  by  the  critical  state  line,  at  the  tip  of  the  advancing  cut.  The 
equation  forms  of  these  two  criteria  will  be  discussed  later. 


Finite  Element  Analysis 

The  analysis  will  typically  be  nonlinear  due  to  a)  the  nonlinear  nature  of  the 
constitutive  equations  that  are  used  to  model  soil  behavior  (for  the  purpose 
of  this  study  a  linear  law  is  used,  however,  in  general  soil  behavior  is 
nonlinear),  b)  the  cutting  criterion  for  the  soil  at  the  tip  of  the  plow,  and  c)  the 
frictional  behavior  occurring  at  the  soil-plow  interface.  A  modified  Newton- 
Raphson  solution  scheme  will  be  used  to  handle  the  nonlinearities.  It  is 
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called  a  "modified"  scheme  because  approximations  are  introduced  into  the 
Jacobian.  These  approximations  have  no  effect  on  accuracy  (as  long  as 
convergence  is  achieved)  but  do  slow  down  the  rate  of  convergence.  The 
desirability  of  introducing  these  approximations  is  explained  as  the  analysis  is 
developed. 

The  iterative  estimate  for  the  vector  of  node  point  degrees  of  freedom 
(displacements  in  the  conventional  displacement  approach  adopted  here  for 
the  undrained  problem)  at  iteration  I  is  written  as  {u}d).  The  Newton- 
Raphson  correction  to  the  I-l  estimate  that  gives  the  Rf*  estimate  is  (the  initial 
estimate  {u}(0)  is  taken  to  be  zero). 


{u)(I)  =  {u}a-l)-[j]a-l)-^R)(I-l)  (8) 

Where  is  the  residual  vector  for  the  finite  element  equations  as  given 

by  the  (I-l)  solution.  The  approximate  Jacobian  is  denoted  by  the  matrix  [J]. 
The  elements  of  the  Jacobian  are  the  first  derivatives  of  the  residuals  with 
respect  to  the  node  point  displacements;  they  are  evaluated  using  the  (I-l) 
estimate  for  the  solution. 

Whether  or  not  convergence  will  occur  and  if  so  the  speed  of  convergence 
depends  in  part  upon  the  accuracy  of  the  approximate  Jacobian  [Jj.  However, 
if  convergence  does  occur,  obtaining  the  correct  solution  (leaving  aside  any 
questions  of  non-uniqueness)  only  depends  upon  using  the  correct  expression 
for  the  residual  vector. 

The  finite  element  equations  are  obtained  using  Galerkin's  weighted  residual 
method  applied  to  the  negatives  of  the  equations  of  motion,  eq(4)  (the 
negative  sign  is  used  to  preserve  the  sign  convention  used  in  a  minimum 
potential  energy  formulation).  The  two  weighted  residuals  associated  with 

node  n  are  (the  same  base  functions  <J>n(x,y)  are  used  to  approximate  both 
displacement  components  ux  and  Uy;  the  base  function  On,  for  node  n,  is 
made  up  of  all  the  element  shape  functions  connected  to  node  point  n;  the 
subscripts  2n-l  and  2n  give  the  locations  of  the  two  rows  in  the  global  residual 
vector): 


dx 


XX 


dx 


dXxy 

dy 


-  P  Vo^  ]  On  dx  dy 


+ 


ay 


-p  Vo2 


9^Uy 

^x2 


]  Ondxdy 


(9) 
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In  the  above  equations,  the  terms  containing  derivatives  of  the  stresses  are 
integrated  by  parts.  This  integration  by  parts  produces  element  interface 
terms  to  cancel  the  residuals  in  the  interface  stress  equilibrium  condition 
(continuity  of  the  traction  vector  across  element  boundaries)  which  are 
implicitly  present  in  the  above  residual  expressions,  the  results  are: 

^2n-l  ^  ”  L  T^xx  dy  P  8x2  J  dx  dy 

(10) 

’^2n=  'W  '^'‘yy  +pvo2  -gjf -I'nJdxdy 

In  a  conventional  finite  element  analysis,  the  integrations  in  the  above 
expressions  are  carried  out  element  by  element  and  give  rise  to  the  element 
contributions  to  the  global  residual  vector  (these  element  contributions  are 
often  called  the  element  load  vectors).  Approximate  (reduced)  integration  is 
introduced  into  this  process.  For  a  four  node,  bi-linear,  isoparametric 
element  the  stress  terms  are  integrated  using  one  point  integration  at  the 
element  center  (i.e.,  at  the  origin  of  the  isoparametric  coordinates  for  the 
element).  This  one  point  integration  is  used  to  avoid  element  "lock-up"  due 
to  the  nearly  incompressible  behavior  of  undrained  soil  problems.  Because 
this  reduced  integration  may  introduce  hour-glassing  into  the  solution,  hour¬ 
glass  control  will  be  used  (e.g.,  see  [4]).  The  inertia  effects  are  integrated  using 
quadrature  points  only  at  the  nodes;  this  process  leads  to  a  lumped  mass 
idealization  which  has  proven  advantageous  in  many  problems  in  dynamics. 

This  partitioning  of  the  integration  into  two  parts,  leads  to  the  assembly  of  the 
finite  element  equations  by  separately  considering  element  contributions  and 
node  contributions.  The  contribution  of  an  element  will  be  contained  in  the 
element  matrices  (the  element  residual  and  Jacobian  matrices,  alternatively 
denoted  as  the  element  load  and  stiffness  matrices). 

From  the  stress  terms  of  the  above  equations  the  components  of  the  element 
residual  vector.  Re,  are  found  (where  the  index  j  runs  from  1  to  4,  c  denotes 
the  element  center,  Nj  are  the  bi-linear  shape  functions  for  the  isoparametric 
element,  A  is  the  area  of  the  element,  and  "Ic"  denotes  evaluation  at  the 
element  center): 


f  8Nj  ^Nj-i 

^2j-l  =  ^  I  +  ^xy  0y 

r  8Nj  3Nj  1 

Re2j  =  At  Xxy  +  tyy  ^ 


(11) 
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It  is  convenient  to  express  the  above  equations  in  matrix  form.  Denote  the 
derivatives  of  the  shape  function  Ni  with  respect  to  x  and  y  as  Fi  and  Gi 
respectively.  The  vector  of  element  residuals  (<Re>  =  <Rei,  Re2,  Res,  ..., 
Re8>)  can  then  be  written  as: 


{Re}  =  A  [B]cT{a]c 


(12) 


Where  (g)  is  the  stress  components  Tij  written  in  vector  form  and  the  matrix 
[B]  is: 


Fi  0  F2  0  F3  0  0 


The  strain  vector  {£}  (<£>  =  <£xx/  Eyy/  Txy>)  can  be  expressed  in  terms  of  the  [B] 
matrix  and  the  vector  of  node  point  displacements  (<U>  =  <Ux|,  Uyp  Ux2/ 
. uy4>): 


{£}  =  [B]  {U} 


(14) 


For  infinitesimal  strains,  the  increment  in  strain  resulting  from  a 
displacement  increment  is  given  by  an  equation  similar  to  eq  (14). 

The  element  Jacobian  (stiffness  matrix  for  a  linear  problem)  is  found  by 
differentiating  the  residual  vector,  eq  (12),  with  respect  to  the  node  point 
displacements.  Before  this  can  be  done  the  stresses  must  be  explicitly 
expressed  as  a  function  of  the  strains.  This  process,  for  a  steady-state  problem, 
presents  some  difficulty.  The  problem  is  that  the  history  dependence  on  the 
strain,  for  a  steady  state  problem,  has  been  converted  to  a  space  dependence, 
see  eq(7).  The  result  is  that  the  stress  at  a  point  such  as  Q  in  Figure  3  is  a 

function  of  the  strains  at  points  Qi,  Q2,  Q3 . Q.  Because  of  this 

dependency,  the  residual  for  element  Q  is  a  function  of  not  only  the 
displacements  of  the  nodes  connected  to  Q,  but  also  of  all  the  displacements  of 
the  nodes  connected  to  the  elements  Qi  to  the  right  of  Q.  If  the  nodes  are 
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numbered  in  the  usual  manner  for  a  conventional  finite  element  analysis  so 
as  to  minimize  the  bandwidth  of  the  equations,  this  dependency  may  result 
in  the  exact  Jacobian  not  being  banded.  In  order  to  avoid  this  possibility,  an 
approximate  Jacobian  is  used.  (Further  study  may  yield  a  more  effective  way 
of  handling  this  problem.)  The  approximate  Jacobian  is  developed  by  writing 
an  approximation  to  eq  (7)  in  the  form: 


{a(Q)}  -  {a(Q)}<l-l)  +  [D(Q)](M)  {Ae(Q)}  (15) 

In  the  above  equation  the  stress  at  the  center  of  a  given  element  (such  as 
point  Q  in  Figure  3)  is  approximated  in  terms  of  the  value  predicted  using  the 
strains  from  the  previous  iteration  (I-l),  for  points  Qi,  Q2,  Qs . Q/  plus  a 


change  induced  by  a  change  (during  iteration  I)  in  strain  at  point  Q.  The 
incremental  change  in  stress  resulting  from  the  incremental  change  in  strain 
at  point  Q  is  predicted  using  the  tangent  stiffness  matrix  [DJd-l)  (i.e.,  {Aa}  =  [D] 
{A8}).  The  tangent  stiffness  matrix  is  given  by  the  algorithm  that  evaluates 
the  material  model  used  to  represent  the  inelastic  behavior  of  the  soil  (note 
that  for  general  inelasticity  the  matrix  [D]  will  not  be  symmetric).  The 
expression  in  eq  (15)  is  an  approximation  as  it  neglects  the  change  in  stress  at 

point  Q  induced  by  the  changes  in  the  strains  at  points  Qi,  Q2,  Qs .  As 

noted  previously  as  long  as  convergence  is  achieved  the  use  of  an 
approximation  Jacobian  will  have  no  affect  on  the  accuracy  of  the  converged 
solution.  The  net  effect  of  this  step  is  to  approximate  the  space  dependence  in 
eq  (7)  by  "successive  substitution"  and  not  by  a  true  Newton-Raphson 
procedure  (for  some  of  the  other  nonlinear  aspects  of  the  problem  the  correct 
contributions  to  the  Jacobian  are  used,  so  the  overall  analysis  will  be  a 
mixture  of  successive  substitution  and  true  Newton-Raphson). 

Introducing  eqs  (14)  and  (15)  into  eq  (12)  and  differentiating  with  respect  to  the 
node  point  displacements  leads  to  the  usual  tangent  stiffness  matrix 
(approximation  to  the  true  Jacobian): 

Qe]=A[B]cT[D]c[B]c  (16) 

The  above  discussion  concerning  the  use  of  the  sequence  of  points  Qi,  Q2, 

Q3 . Q  to  evaluate  history  dependency,  suggests  that  the  use  of  a  mesh 

where  the  grid  lines  in  the  x  direction  are  straight  and  parallel  to  the  x  axis  is 
advantageous  (see  Figure  4  for  an  example  of  a  very  course  grid).  With  this 
property  the  points  Qi  can  all  be  located  at  element  centers  (if  such  a  mesh  is 
not  used  an  interpolation  scheme  must  be  used  to  find  the  strains  at  points 
Ql/  Q2/  Qs-  /  which  all  lie  at  the  same  distance  y  from  the  axis  of  symmetry  as 
does  point  Q,  in  terms  of  the  strains  at  the  centers  of  the  surrounding 
elements).  For  simplicity  the  remainder  of  the  analysis  is  restricted  to 
meshes  of  the  type  shown  in  Figure  4. 
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a)  portion  of ‘deformed  grid 


■< - type  c - 1 - >- 


Figure  4.  Deformed  and  undeformed  grids  for  2-D  plow  problem 


Considering  element  "n"  of  Figure  4,  denote  the  x,y  coordinates  of  the 
connecting  nodes  as  Xi,yi  (with  node  ni  selected  to  be  the  upper  right  one)  ; 

the  isoparametric  coordinates  are  and  T|i=(l,l,  -1,-1).  The  well 

known  isoparametric  transformation  can  be  easily  specialized  for  the  special 
element  type  of  Figure  4  and  showm  to  yield: 

Ni,=  l/4 

Fi,=^i/a2  (17) 

ai/a2)/2Ay 

Where  (the  element  area  is  A): 
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ai  =  XJ+X2-X3-X4 


a2  =  X1-X2-X3+X4 


as  =  X1-X2+X3-X4 


(18) 


Ay=yi-y4 

A  =  Ay  as  /2 

For  later  use  the  values  of  Fi,  Gi  and  Ni  at  the  corner  k  of  the  element  are 
needed  (5ik  is  the  Kronecker  delta,  i.e.,  the  identity  matrix): 


Nik  =  5ik 


Fik  =^i  (1  +  TliTlk)/(a2  +  Tlk  as )  (19) 

Gik  =[r[i  (1  +  ^i^k)  -  (1  +  TliTlk)  (ai  +  Tik  as  )/(a2  +  rik  as  )]/2  Ay 

The  element  Jacobian  and  residual  matrices  (eqs  12  and  16)  are  assembled  in 
the  usual  way  using  the  "direct  stiffness"  concept;  to  this  must  be  added  the 
contributions  from  the  integration  of  the  inertia  terms  in  eq  (10)  using 
quadrature  points  placed  at  the  nodes. 

Consider  the  node  "m"  in  Figure  4  surrounded  by  the  area  with  nodes  "a" 
and  "b"  directly  to  the  left  and  to  the  right.  A  central  finite  difference  operator 
is  used  to  approximate  the  second  derivatives  with  respect  to  x  at  point  m: 


a2ux 

=  c.|  Uxa  +Co  Uxm  +C|  Uxb 


a^Uy 

=  C.1  Uya  +Co  Uy^  +C^  Uy^ 

where 

_ 2 

^  (Xm-Xa)  (Xb-Xa) 

-2 

(Xm-Xa)  (Xb-Xm) 


2 

(xb-Xm)  (xb-Xa) 


(20) 


(21) 
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The  residual  contributions  of  the  last  terms  in  eq  (10)  at  node  m,  can  now  be 
evaluated  (ux  and  Uy  are  the  I-l  estimates): 


Rnode2m-l  =  «  (c.^  Ux^  +0^  Ux^  +C|  Ux^) 
Rriode2m  =  W  (c.j  Uy^  +Cq  Uyj^  +C|  Uy^,) 


(22) 


where 


CO  =  Am  p  Vo2 

The  contributions  to  the  Jacobian  (Jnodei^j  denotes  the  contribution  to  the  i,j 
term  in  the  matrix  [J])  are  found  by  differentiating  the  residual  contributions 
with  respect  to  the  displacement  components: 

}node2m-l,  2a-l  =  C_;j 

Jnode2m-l,  2m-l  =  ®  C^ 

Jnode2m-l,  2b-l  =  CO  C|  (23) 


Jnode2m,  2a  =  CO  C.| 

Jnode2m,  2m  =  CO  Cq 
Jnode2m,  2b  =  CO  C| 

Special  consideration  is  required  when  "m"  is  a  node  on  either  the  left  or 
right  boundary  of  the  mesh  (i.e.,  either  no  node  point  "a"  to  the  left  or  "b"  to 
the  right  exist).  It  is  required  that  the  horizontal  dimensions  of  the  soil  mass 
being  analyzed  be  selected  large  enough  that  at  the  very  left  the  soil  has  come 
to  rest  and,  thus,  there  is  no  acceleration  while  at  the  right  the  soil  has  as  yet 
not  felt  the  presence  of  the  advancing  plow  and,  thus,  no  acceleration  exists. 
Hence,  the  node  contributions  to  the  residual  and  Jacobian  matrices  are  zero 
at  these  points. 

After  the  contributions  of  all  elements  and  node  points  to  the  global  Jacobian 
and  residual  matrices  have  been  assembled,  the  resulting  matrices  must  be 
modified  to  account  for  boundary  conditions  before  the  solution  of  eq  (8)  is 
undertaken.  In  Figure  4  the  several  different  boundary  types  are  identified  by 
the  letters  "a"  through  "f".  Type  "a"  is  the  line  of  symmetry  where  Uy=0  (the 
specification  of  displacement  boundary  conditions  is  done  in  the  usual  way 
for  finite  element  analyses).  Boundary  "b"  is  required  to  be  far  enough 
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removed  that  the  soil  has  not  yet  felt  the  disturbance  of  the  advancing  plow, 
hence,  ux=0  and  uy=0.  For  the  plowing  of  a  large  soil  mass,  boundary  "c"  is 
required  to  be  far  enough  removed  that  the  soil  does  not  detect  the  passing 
plow,  thus,  ux=0  and  uy=0.  Boundary  "d"  is  required  to  be  far  enough  to  the 
left  that  the  soil  has  reached  steady  state  behavior,  thus 


Using  eq  (2) 


3ux 

=  0 

at 

aUy 

at 

=  0 

this  ] 

gives 

aux 

=  0 

ax 

auy 

ax 

=  0 

(24) 


(25) 


The  derivative  boundary  conditions,  for  a  point  such  as  "j"  on  boundary  d, 
are  specified  using  two  point  finite  difference  operators: 


Xk-Xj 


Xk-Xj 


=  0 


(26) 


UXj-UXk  -  0 

uyj'Uyk=  0 


(27) 


Because  the  soil  model,  the  interface  condition  along  the  plow's  surface  and 
the  cutting  condition  specification  for  the  point  at  the  tip  of  the  plow,  in 
general,  all  destroy  the  symmetry  of  the  Jacobian  matrix  no  attempt  is  made  to 
preserve  symmetry  in  the  implementation  of  eq  (27).  The  implementation 
of  eq  (27)  is  easily  accomplished  by  replacing  the  finite  element  equilibrium 
equation  2j-l  by  the  equation  that  uxj-uxk=0.  Similarly  equation  2j  is  replaced 
by  the  equation  uyj-uyj^=0.  The  residuals  just  become  the  differences  between 
the  two  displacements  as  estimated  in  the  previous  iteration;  the  Jacobian,  the 
derivative  of  the  residual,  has  entries  of  zeros  and  ones  in  the  appropriate 
columns,  i.e.. 
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(28) 


R,i  =  Uy,(I-«-Uy,a-l) 
and  (all  other  elements  are  zero) 


J2H/2k-l  = 
J2j.2j  =  1 


(29) 


J2j/2k  - 

Along  the  boundary  "e"  it  is  assumed  that  the  trench  produced  by  the  plow 
remains  open  and,  thus,  the  surface  traction  acting  on  the  soil  is  just  the 
water  pressure  if  the  plowing  is  below  the  surface  of  water  and  zero  if  it  is  in 
air;  the  resulting  node  point  forces  are  added  to  the  appropriate  rows  of  the 
residual  in  the  usual  fashion. 


Boundary  segment  f  is  the  interface  (excluding  node  pc)  between  the  soil  and 
the  plow.  Denote  the  equation  of  the  surface  of  the  plow  by  y=g(x’^)  where  x* 
is  measured  from  the  plow  tip,  i.e,  x*=x-uc,  (recall  that  at  the  time  of  interest 
T,  the  location  of  the  coordinate  axis  is  specified  so  that  the  tip  of  the  plow  is 
at  x=Uc,  see  Figures  2 , 3,  4  and  5).  Along  the  soil-plow  interface  the  frictional 
relationship  between  the  effective,  normal  and  tangential  interface  stresses 
must  be  enforced.  This  condition  is  approximated  in  terms  of  the 
"generalized  forces"  at  the  several  nodes  lying  along  the  interface,  such  as 
node  "i"  in  Figure  5.  Both  the  virtual  work  and  the  minimum  potential 
energy  interpretation  of  the  finite  element  equations  can  be  used  to  give  the 
interpretation  of  the  negatives  of  the  node  point  residuals  (-R^j  ^  and  -R^j)  as 

external  reactions  applied  at  the  node  (these  of  course  should  be  zero  at 
interior  and  stress  free  boundary  nodes).  Thus  in  Figure  5; 


Fxi  =  -R,,. 


2i-leff 


Fyi  - 


(30) 
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Figure  5.  Components  of  generalized  force  at  node  i, 
due  to  action  of  soil  on  plow 

The  "eff"  subscript  indicates  that  the  residuals  must  be  those  calculated  using 
effective  soil  stresses  (as  opposed  to  using  total  stress  as  is  the  case  for  all  other 
nodes  where  equilibrium  equations  are  being  generated).  This  mixed  use  of 
total  and  effective  stress  based  residuals  is  accomplished  by  replacing  those 
rows  in  the  element  residual  and  Jacobian  matrices,  which  correspond  to 
nodes  along  the  interface,  by  rows  calculated  using  the  effective  stresses  (i.e., 
in  eq  12  effective  stresses  are  used  and  in  eq  16  the  effective  tangent  stiffness 
matrix  is  used,  that  is  it  is  not  augmented  by  the  effective  bulk  modulus  Keff)- 
Equation  30  is  now  used  to  find  the  normal  (n)  and  tangential  components  (t) 
of  generalized  force: 


(31) 


Assuming  that  the  interaction  of  the  plow  and  the  soil  can  be  described  by 
Coulomb  friction  with  a  coefficient  of  friction  of  "f"  (the  direction  of  Ftj  must 
be  as  shown  in  Figure  5  in  order  to  oppose  the  motion  of  the  plow): 

Fti-fFni=0  (32) 


Introducing  eq  (31)  into  eq  (32)  gives: 

R-.  ,  [sin(e.)  +  f  cos(0.)]  -  R,.  „  [cos(03  -  f  sin(0  )]  =  0  (33) 
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This  equation  replaces  the  (2i-l)th  finite  element  equation.  The  new  residual 
is  the  left  hand  side  of  eq  (33)  (where  R2J  ^  and  R^j  ^^are  the  old  rows  in  the 

residual  vector  calculated  using  the  effective  soil  stresses  not  the  total) 
evaluated  using  the  (I-l)  estimate  of  the  solution.  Differentiating  eq  (33)  with 

respect  to  the  displacement  components  (and  neglecting  the  dependence  of  0. 

on  Uc)  gives  the  new  row  in  the  global  Jacobian  matrix  as  a  combination  of 
the  old  rows  where  the  R's  of  eq  (33)  are  replaced  by  J's.  The  sin  and  cos  of 
the  angle  9.  are  calculated  from  geometry: 

sin(0p 


cos(6p 

Where  g'  is  the  derivative  of  g  with  respect  to  x*.  The  values  of  g  and  g’  are 
evaluated  at  the  deformed  location  of  node  "i",  i.e.,  at  (x*=xi  +  Uxj-  uc),  where 
the  values  of  ux,  and  uc  are  estimated  using  the  results  of  the  previous 
iteration.  Because  this  dependency  of  6.  on  the  solution  is  ignored  when 

forming  the  new  row  for  the  Jacobian  it  is  an  approximation  to  the  true 
Jacobian  (as  noted  previously  such  an  approximation  has  no  affect  on  the 
accuracy  of  the  converged  solution). 

In  addition  to  the  frictional  law,  along  the  soil-plow  interface,  the  conformity 
of  the  soil  to  the  plow  profile  must  be  enforced  (it  is  assumed  that  no 
separation  occurs).  Consider  node  "i"  on  the  plow-soil  interface,  the 
compatibility  condition  between  the  plow  and  soil  requires  that: 

Uyj  -  g(Xi+Uxf Uc)  =  0  (35) 

This  equation  replaces  the  (2  i)l^^  finite  element  equation.  The  residual  is  the 
left  hand  side  of  eq  (35)  (where  the  estimates  from  iteration  I-l  are  used  for 
uxj  and  uc).  The  new  row  for  the  Jacobian  is  found  by  differentiating  the 
residual  with  respect  to  the  displacement  components  (again  the  dependence 
on  Uc  is  ignored  in  this  process): 


.  .  3' _ 

Vl+(g')2 

(34) 

1 


J 


=  -g 


(36) 
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The  above  process  assures  that  the  deformed  shape  of  the  soil  has  the  precise 
contour  of  the  plow  (assuming  that  the  plow  can  be  modeled  as  rigid).  The 
precise  fitting  of  the  deformed  soil  to  the  plow  profile  may  be  judged  to  be  a 
large  deformation  effect  which  could  possibly  be  approximated  by  neglecting 
the  difference  between  x  and  x*.  However,  because  the  precise  specification  is 
quite  a  simple  matter  it  is  included  in  the  analysis. 

At  node  pc  uy=0.  In  addition,  the  "cutting  criterion"  for  the  soil  must  be 
specified.  The  determination  of  what  the  cutting  state  should  be  is  a  problem 
in  fracture  mechanics  which  to  the  Author's  knowledge  has  not  been 
addressed  for  soils.  In  this  work  two  simple  assumptions  for  the  cutting  state 
are  discussed. 

The  first  is  that  the  soil  being  cut  has  reached  a  limiting  value  of  effective 
tensile  stress,  Ofaii,  across  the  surface  being  cut.  This  assumption  is 
implemented  by  requiring  that  the  effective  normal  (to  the  cut)  stress  Xyygff 

should  be  equal  to  Ofaii  at  point  pc-  This  specification  replaces  the  finite 
element  equilibrium  equation  for  the  x  direction  at  node  pc  (the  equilibrium 
equation  contains  the  generalized  force  exerted  on  the  soil  by  the  tip  of  the 
plow,  this  force  is  not  known  a  priori): 

^yyeffpc  -  =  0  (37) 

The  quantity  Xyygffp^  is  the  effective  normal  stress  at  point  pc  as  predicted 

using  the  finite  element  approximation  in  element  c.  The  prediction  of  the 
state  of  stress  at  point  pc  in  element  c  requires  that  the  history  dependency  be 
evaluated  using  the  strain  states  for  all  the  corresponding  points  in  the  row  of 
elements  to  the  right  of  c.  This  process  was  previously  discussed  for  the 
prediction  of  the  states  of  stress  at  the  element  centers. 

The  (2pc-l)  residual  is  the  left-hand  side  of  eq  (37)  as  predicted  using  the 
solution  from  the  (I-l)  iteration.  The  corresponding  row  in  the  Jacobian  is 
found  by  differentiating  the  residual  with  respect  to  the  displacement 
components  (use  is  made  of  eq  15): 

hpc-U  "  ^^'^eff  ^n,j  (38) 

Where  Dijgff  are  the  tangent  stiffness  properties  for  the  effective  stress  at  the 
location  of  the  third  node  (i.e.,  pc)  of  element  c.  The  components  of  the  B 
matrix  (eq  13)  are  found  using  eq  (19)  with  k=3;  the  repeated  index  n  is 
summed  from  1  to  3,  and  the  subscript  j  in  the  Jacobian  refers  to  the  eight 
degrees  of  freedom  for  element  c  (i.e.,  2mi-l  and  2mi,  where  mi  are  the  four 
nodes  for  element  c). 
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The  second  possibility  for  a  cutting  criterion  at  node  pc  was  suggested  by 
Professor  Kutter.  It  assumes  that  the  soil  is  in  a  state  of  failure  at  point  pc  as 
described  by  critical  state  soil  mechanics.  This  condition  assumes  that  cutting 
occurs  when  the  effective  stress  state  for  the  soil  lies  on  the  critical  state  line 
in  1,  J  space  (I  is  three  times  the  mean  effective  stress  and  J  is  the  second  stress 
invariant  for  the  deviatoric  stress),  i.e., 

J/I-M/V27  =  0  (39) 

Where  M  is  the  slope  of  the  critical  state  line  in  p,  q  space.  The 
implementation  of  this  cutting  condition  replaces  the  2pc-l  residual  with  the 
left-hand  side  of  eq  (39)  evaluated  at  the  third  node  of  element  c.  The  new 
row  in  the  Jacobian  is  found  by  differentiating  the  residual  (see  the  discussion 
of  eq  39  for  the  meaning  of  J,  etc.): 


V3 


where 

2 

Pin  =  X 

i=l 


5  I 

i=l 


(40) 


(41) 


With  the  boundary  condition  modifications  of  the  finite  element  equations 
completed,  the  equations  are  ready  for  solution  to  give  the  Newton-Raphson 
correction  to  the  node  point  displacements  for  iteration  I,  see  eq  (8).  As  noted 
previously  the  simultaneous  equations  are  banded  but  non-symmetric. 
Iteration  must  be  continued  until  convergence  is  achieved.  Convergence  is 
determined  by  the  convergence  of  the  plow  force  (described  below). 

It  is  of  primary  interest  to  determine  the  force  required  to  propel  the  plow 
through  the  soil.  This  force  is  calculated  as  2  times  (to  account  for  the  fact  that 
symmetry  is  used  in  the  modeling  of  the  plow)  the  sum  of  all  the  x 
components  of  the  total  generalized  forces  acting  at  the  interface  nodes  along 
the  soil-plow  interface  (including  point  pc)  minus  any  force  due  to  water 
pressure  acting  on  the  back  of  the  plow. 


Example 
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The  theory  presented  in  the  previous  section  has  been  implemented  in  a 
simple,  non-production,  finite  element  code  intended  for  preliminary 
evaluation  of  the  one-step,  steady-state  analysis  concept.  The  analysis  was 
implemented  for  both  plane  stress  and  plane  strain  conditions  for  a  two- 
dimensional  idealization  of  the  plow  problem.  No  attempt  was  made  to 
accurately  model  the  soil’s  inelasticity  nor  the  pore-water  pressure 
redistribution.  Instead  a  linear  viscoelastic  model  was  used  to  model  the  soil, 
pore-water  system  (this  choice  was  made  in  order  to  simplify  and,  thus, 
expedite  this  preliminary  implementation).  A  linear  viscoelastic  model  can 
qualitatively  model  the  phenomena  of  soil  inelasticity  and  pore-water 
pressure  redistribution  but  can  not  do  so  quantitatively.  Because  it  can 
qualitatively  model  these  phenomena  it  is  felt  that  its  inclusion  in  the  trial 
implementation  will  demonstrate  whether  or  not  the  one-step,  steady-state, 
inelastic  analysis  method  is  capable  of  handling  these  phenomena. 

Because  of  the  gross  idealizations  of  two-dimensional  geometry  and  linear 
viscoelastic  behavior  of  the  soil,  pore-water  system,  no  quantitative 
significance  should  be  given  to  the  results  that  are  presented  below.  The 
example,  however,  will  demonstrate  the  ability  of  the  analysis  method  to 
handle  either  steady-state  dynamic  or  quasi-static  conditions;  time-dependent, 
inelastic  material  behavior;  interface  friction  between  plow  and  soil; 
prescription  of  a  failure  condition  for  the  soil  at  the  tip  of  the  plow;  and  the 
ability  to  accuratdy  model  arbitrarily  shaped  plows. 

Two  possible  "cutting  conditions"  were  suggested  in  the  previous  section. 
For  the  purpose  of  this  example  the  imposed  cutting  condition  was  that  the 
normal  stress  in  the  soil  just  ahead  of  the  plow  tip  is  zero  (assumed  zero 
tensile  strength  of  the  soil).  While  this  condition  may  be  a  gross 
simplification  of  the  actual  situation,  it  is  thought  that  comparisons  of 
competing  plow  designs  would  still  be  meaningful  and  that  quantitative 
predictions  of  plow  force  would  still  be  relatively  accurate  (if  the  restrictions 
of  two-dimensional  geometry  and  idealized  linear  viscoelastic  soil  behavior 
were  removed). 

As  an  illustration  of  the  potential  capabilities  of  the  one-step,  steady-state 
analysis  concept,  the  following  problem  was  analyzed.  A  plane  stress  analysis 
was  performed  of  a  plow  moving  at  constant  velocity  Vq  through  a  layer  of 
soil,  see  Figure  6.  For  the  purpose  of  simplicity  of  data  generation,  the  shape 
of  the  plow  was  taken  to  be  a  simple  analytical  curve  of  ybiade  =  t/2 
where  -L  <  x*  <  0  (the  coordinate  x*  is  measured  from  the  tip  of  the  plow,  see 
discussion  in  the  previous  section).  The  length  of  the  blade  is  denoted  by  L, 

the  blade  thickness  is  t  and  a  is  a  parameter  that  controls  blade  shape.  As  was 
stated  it  was  for  the  sake  of  convenience  that  a  simple  equation  was  assumed 
for  the  shape  of  the  blade;  it  is  in  fact  a  very  simple  matter  to  accommodate 
arbitrarily  shaped  plows. 


21 


The  assumed  linear  viscoelastic  properties  of  the  soil,  pore-water  system  are 
described  by  two  relaxation  functions,  one  in  shear  and  the  other  in  volume 
change,  i.e., 

OG(t)  =  Goo  +  (Go  -  Goo)  e-«Gt 
OB(t)  =  Boo  +  (Bo  -  Boo)  e-“Bt 

The  initial  shear  and  bulk  moduli  are  subscripted  by  "o",  while  the  values  at 
infinite  time  are  subscripted  by  The  rates  of  relaxation  are  controlled  by  ac 
and  ttB- 

Using  symmetry,  the  analysis  can  be  restricted  to  half  of  the  soil  mass,  see  the 
lower  part  of  Figure  6.  The  boundary  condition  at  the  outer  edge  of  the  soil 
mass  (the  tank  wall  for  the  Navy  test,  see  [5,6])  was  taken  to  be  a  fixed 
boundary. 

Several  analyses  were  run  in  order  to  determine  the  placement  of  the  right 
boundary  so  that  it  was  sufficiently  far  upstream  of  the  plow  not  to  feel  any 
effect  of  the  advancing  plow.  The  boundary  condition  at  this  edge  was  taken 
to  be  fixed.  Several  analyses  were  run  in  order  to  determine  the  placement  of 
the  left  boundary  so  that  it  was  sufficiently  far  downstream  of  the  plow  to 
have  reached  a  uniform  state.  The  condition  of  uniformity  was  specified  by 
setting  the  derivatives  of  the  displacement  components  with  respect  to  x  to  be 
zero  at  each  of  the  node  points  along  the  left  side  boundary. 

The  boundary  condition  behind  the  plow  and  along  the  center-line  of  the  soil 
is  one  of  a  stress  free  surface.  The  boundary  condition  ahead  of  the  plow  is 
the  condition  of  symmetry.  At  the  interface  of  the  plow  and  the  soil  (except  at 
the  plow  tip)  the  conditions  are  that  the  soil  conforms  to  the  shape  of  the 
plow  and  the  generalized  nodal  shear  forces  equal  the  coefficient  of  friction  f 
times  the  generalized  nodal  normal  forces. 
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a)  Problem  Configuration 


a.'.'.'.' 

. . . 

- 

b)  Configuration  for  Analysis 


Figure  6.  Simple  plow  problem 

Finally  at  the  very  tip  of  the  plow,  the  displacement  Uy  is  required  to  be  zero 
and  the  cutting  condition  is  specified. 

The  total  force  Fq,  required  to  move  the  plow  at  constant  velocity  through  the 
soil,  is  calculated  by  summing  the  x  components  of  the  node  point 
generalized  forces  for  the  nodes  along  the  plow-soil  interface  (including  the 
tip  of  the  plow). 

The  coupling  introduced  into  the  simultaneous  equations  (for  the  steady- 
state,  history  dependent,  finite  element  analysis)  by  the  interchange  of  time 
and  space  when  evaluating  the  history  dependence  of  the  soil  properties  was 
handled  by  moving  the  coupling  terms  to  the  right-hand  side  of  the  equations 
and  approximating  them  by  iteration.  (Thus,  the  banded  nature  of  the  finite 
element  equations  was  preserved.)  Iteration  was  of  course  already  necessary 
to  model  the  nonlinearities  of  the  problem  (material  inelasticity,  the 
frictional  interface  condition  and  the  enforcement  of  the  cutting  condition  at 
the  tip  of  the  plow).  As  was  previously  noted,  the  equations  are  non- 
symmetric  (this  is  also  usually  true  for  a  time  marching  analysis  if  the  correct 
Jacobian  is  used  for  the  inelastic  material). 
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Most  of  the  parameters  describing  the  example  problem  were  selected  with 
the  Navy's  laboratory  test  [5,6]  in  mind,  however,  because  the  assumed 
viscoelastic  behavior  for  the  soil  can  at  best  qualitatively  represent  the 
inelasticity  of  the  soil  and  the  movement  of  pore  water,  no  effort  was  made  to 
correlate  this  aspect  of  the  problem  description  to  the  actual  test  conditions. 
The  selected  parameters  were:  viscoelastic  soil  properties  of  Goo=50  psi,  Go=500 
psi,  ocg=30  sec'^  and  Boo=Bo=335,000  psi;  soil  density  of  .081  Ib/in^;  blade  shape 

and  dimensions  of  L=6.5  in,  t=3/8  in  and  a=0.7;  the  dimensions  of  the  soil 
mass  were  45  inches  to  the  left  of  the  plow,  45  inches  to  the  right  of  the  plow 
and  15  inches  for  the  half  width: 

A  basic  4  element  grid  was  selected  and  successively  refined  (maximum  of 
approximately  2500  elements)  until  the  predicted  value  of  Fq  converged 
(results  from  the  1300  element  and  2500  element  grids  were  nearly  identical). 
From  5  to  10  iterations  were  required  for  convergence  of  the  one-step,  steady- 
state  analysis.  Plow  velocities  of  1,  3, 5  and  7  fps  were  considered;  results  for  5 
fps  are  given  in  Figures  7-10. 

Figure  7  shows  the  deformed  mesh  (displacements  and  plow  thickness 
magnified  by  a  factor  of  5);  all  elements  were  rectangles  in  the  undeformed 
mesh.  Contour  plots  of  the  soil  strains  Jxy  and  Ex  are  given  in  Figures  8  and  9, 
a  contour  plot  of  the  soil  stress  Txy  is  given  in  Figure  10.  These  plots  clearly 
show  that  the  length  dimension  has  been  taken  large  enough  so  as  to  reach 
undisturbed  soil  on  the  right  and  uniformity  conditions  on  the  left.  Given  in 
Figure  11  are  plots  of  the  plow  force  (per  unit  thickness  of  the  plane  stress 
body)  as  a  function  of  plow  velocity.  The  circles  are  for  the  plow 
configuration  of  a=0.7  and  the  case  where  acceleration  of  the  soil  is  included; 
when  the  inertia  is  neglected  the  curve  with  square  symbols  was  obtained. 
The  curve  with  diamond  symbols  is  for  the  case  of  a  differently  shaped  plow 

(a=1.2);  the  two  plow  shapes  are  illustrated  in  Figure  12  (the  vertical  scale  is 
exaggerated  by  about  a  factor  of  7;  the  waviness  of  the  lines  are  due  to  the 
plotting  program). 
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mesh  limits: 
xmin  -5.1500E+01 
xmax  4.5000E+01 
ymin  O.OOOOE+00 
ymax  1.5000E+01 


status  data: 
nodes  used  1440 
elements  1357 


Figure  7:  Deformed  mesh  for  example  plow  problem,  vo=5  fps 


strain  gamma -xy 


min  -1.838E-01 
max  4.2e2E-03 


■  ♦2  .  OOOOE-02 
+1.2000E-02 
♦  4  .  OOOOE-03 
-4.0000E-03 
L,'  -1.2000E-02 

-2.0000E-02 
-2 . 8000E-02 


-3 . 6000E-02 
-4 .4000E-02 
-5.2000E-02 
-6 . OOOOE-02 
-6 . 8000E-02 
-7 .6000E-02 
-8.4000E-02 
-9 .2000E-02 
-1. OOOOE-01 
-1.0800E-01 
-i. 1600E-01 
-1. 2400E-01 
-1.3200E-01 
-1.4000E-01 
-1.4800E-01 
-1.5600E-01 
-1.6400E-01 
-1.7200E-01 
-1.8000E-01 
-1.8800E-01 


Figure  8:  Contour  plot  of  shear  strain  (yxy)  hi  soil,  vo=5  fps 
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strain  epsilon-x 


min  -1.189E-01 
max  2.997E-02 


■I  +3 .3000E-02 
I  +2 .7000E- 02 
I  +2 . lOOOE-02 
,  +1. 5000E-02 

+9  .OOOOE-03 
+  3  .  OOOOE-03 
-3.0000E-03 


-9  .  OOOOE-03 
-1. 5000E-02 
-2 . lOOOE-02 
-2 .7000E- 02 
•3.3000E-02 
-3.9000E-02 
-4 .5000E-02 
-5 . lOOOE-02 
-5 .7000E-02 
-6.3000E-02 
-6.9000E-02 
-7 . 5000E-02 
-8.1000E-02 
-8.7000E-02 
-9 . 3000E- 02 
-9 . 9000E-02 
-1. 0500E-01 
-l.llOOE-01 
-1.1700E-0i 
-1.2300E-01 


Figure  9:  Contour  plot  of  normal  strain  (Ex)  in  soil,  vo=5  fps 


27 


etress  tau-xy 


min  -8,403E+01 
max  3.036E+00 


■+1.4000E+01 
+1. OOOOE+01 
+6 . OOOOE+00 
+2 . OOOOE+00 
-2.0000E  +  00 
-6.0000E  +  00 
-l.OOOOE+01 


-  1 . 4000E  +  01 
-1. 9000E+01 
-2 . 2000E+01 
-2 . 6000E  +  01 
-3  .  OOOOE  +  01 
-3 . 4000E  +  01 
-3 . 8000E+01 
-4 . 2000E+01 
-4 . 6000E+01 
-5 . OOOOE+01 
-5.4000E+01 
-5.8000E+01 
-6 . 2000E+01 

-  6 . 6000E  +  01 
-7 . OOOOE+Oi 
-7 .4000E+01 
-7 . 8000E+01 
-8 . 2000E+01 
-8. 6000E+01 
-9 . OOOOE+01 


Figure  10:  Contoxir  plot  of  shear  stress  (Txy)  in  soil,  vo=5  fps 
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Force 


0) 


distance  along  plow-x* 


Figure  12.  Shapes  of  plow  for  example  problem 

It  is  to  be  emphasized  that  no  attempt  was  made  to  calibrate  the  viscoelastic 
material  model  in  order  to  capture  the  actual  soil  inelasticity  and  pore-water 
pressure  redistribution  nor  to  pick  its  parameters  so  that  Fq  would  agree  with 
experimental  observations  (it  appears  that  this  would  be  possible  but  it  would 
have  very  little  significance). 

What  is  important  to  observe  from  these  results  is  the  predicted  dependence 
of  the  plow  force  Fq  on  the  velocity,  Vq,  the  rather  dramatic  dependence  of  the 
force  upon  the  plow  shape  and  the  importance  of  including  inertia  effects. 
However,  of  even  greater  significance  is  the  observation  of  the  inexpensive 
nature  of  the  analysis.  As  was  mentioned  only  5  to  10  iterations  were 
required  for  convergence.  Of  course,  because  it  is  a  steady-state  analysis,  no 
time  stepping  is  required. 


Conclusions 

From  this  preliminary  study  several  conclusions  can  be  drawn.  Because  of 
the  very  few  iterations  (and  no  time  stepping)  required,  it  should  be  entirely 
feasible  to  perform  a  three-dimensional,  steady-state,  inelasticity  study  of  the 
plow  problem  including  parameter  studies  for  different  plow  shapes  and 
plowing  velocities.  Thus,  it  appears  that  the  one-step,  steady-state  analysis 
procedure  offers  a  viable  alternative  to  a  multi-step,  transient  analysis  of  the 
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plow  problem.  However,  because  no  three-dimensional  code  currently  exists 
for  the  one-step,  steady-state  procedure,  it  would  seem  that  the  best  way  to 
proceed  for  the  present  study  is  to  use  an  available  commercial  code  to 
perform  a  multi-step,  transient  analysis. 

Before  any  attempt  is  made  to  produce  a  production  code  for  the  one-step, 
steady-state  procedure  several  items  require  further  investigation.  The 
failure  (or  cutting)  condition  in  the  soil  at  the  tip  of  the  plow  (i.e.  the  state  of 
the  soil  as  it  is  cut  by  the  plow  tip  to  form  a  new  surface)  is  not  well 
understood  and  is  an  area  that  should  receive  further  theoretical  and 
experimental  study  (this  same  information  is  also  required  for  a  rigorous 
multi-step,  transient  analysis  of  the  plow  problem).  The  affect  of  using  a 
realistic  plasticity  model  for  the  soil  on  the  rate  of  convergence  of  the 
iteration  process  must  be  studied.  Finally,  means  for  the  incorporation  of  the 
flow  of  pore-water  into  the  analysis  must  be  developed. 


References 

1.  Baligh,  M.  B.  (1985),  "Strain  Path  Method",  Journal  of  Geotechnical 
Engineering,  Vol.  Ill,  No.  9,  pp.  1108-1136 

2.  Baligh,  M.  B.,  Azzouz,  A.  S.  and  Chin,  C.T.  (1987),  "Disturbances  Due  to 
'Ideal'  Tube  Sampling",  Journal  of  Geotechnical  Engineering,  Vol.  113, 
No.  7,  pp.  739-757 

3.  Christensen,  R.  (1982),  Theory  of  Viscoelasticity:  An  Introduction,  2nd 
edition.  Academic  Press 

4.  Hughes,  T.  J.  R.  (1987),  The  Finite  Element  Method,  Prentice-Hall  Inc. 

5.  Cable,  S.B.,  Carlisle,  H.,  Skyers,  B.,  and  Taylor,  R.  (1993),  "Hydrostatic 
Pressure  Effects  on  Thin  Blade  Plow  Resistance  in  Dense,  Saturated, 
Cohesionless  Soil",  Report  No.  TM-42-93-06,  Naval  Civil  Engineering 
Laboratory,  Port  Hueneme,  CA. 

6.  Girard,  J.,  and  Taylor,  R.  (1994),  "Blade  Geometry  and  Soil  Permeability 
Effects  on  thin  Blade  Plow  Resistance  in  Dense,  Saturated,  Cohesionless 
Soils  -  Phase  II  Report".  Report  No.  TM-2026-OCN,  Naval  Facilities 
Engineering  Service  Center,  Port  Hueneme,  CA. 


31 


DISTRIBUTION  LIST 


AC  ENGRG  INC  /  BERRY,  WEST  LAFAYETTE,  IN 

ADINA  ENGRG,  INC  /  WALCZAK,  WATERTOWN,  MA 

AEWES  /  UB,  VICKSBURG,  MS 

AEWES  /  PETERS,  VICKSBURG,  MS 

AFWL/NTE  /  BALADI,  KIRTLAND  APB,  TX 

ANATECH  APPLICATIONS  /  CASTRO,  SAN  DIEGO,  CA 

ANATECH  RESEARCH  CORP  /  RASHID,  SAN  DIEGO,  CA 

APPLIED  PHYSICS  TECH  /  SWANSON,  MCLEAN,  VA 

APPLIED  RSCH  ASSOC,  INC  /  HIGGINS,  ALBUQUERQUE,  NM 

ARMY  EWES  /  WES  (NORMAN),  VICKSBURG,  MS 

ARMY  EWES  /  WESIM-C  (N.  RADHAKRISHNAN),  VICKSBURG,  MS 

ASSOCIATED  SCIENTISTS  /  MCCOY,  WOODS  HOLE,  MA 

BING  C  YEN  /  IRVINE,  CA 

BRITISH  EMBASSY  /  ELLIS,  WASHINGTON,  DC 

BYU  /  ROLLINS,  PROVO,  UT 

CAUF  INST  OF  TECH  /  PASADENA,  CA 

CALTRANS  OFFICE  OF  RESEARCH  /  HOLLAND,  SACRAMENTO,  CA 
CATHOLIC  UNTV  /  CE  DEPT  (KIM)  WASHINGTON,  DC 
CENTRIC  ENGINEERING  SYSTEMS  INC  /  TAYLOR,  PALO  ALTO,  CA 
CHALMERS  UNIVERSITY  OF  TECHNOLOGY  /  TEPFERS,  412  74  GOTEBORG, 
CHEUNG  AND  ASSOCIATES  /  CHEUNG,  IRVINE,  CA 
COLORADO  SCHOOL  OF  MINES  /  GOLDEN,  CO 
COLORADO  ST  UNIV  /  FORT  COLLINS,  CO 
COMPUTATIONAL  MECHANICS  /  BREBBIA  SOUTHAMPTON, 

CORNELL  UNIV  /  ITHACA,  NY 

COUNTY  OF  VENTURA  /  TAKAHASHI,  VENTURA,  CA 
CRREL  /  KOVACS,  HANOVER,  NH 
CSU  CHICO  /  ARTHUR,  CHICO,  CA 
CSU  FULLERTON  /  RAMSAMOOJ,  FULLERTON,  CA 
DAMES  &  MOORE  /  LOS  ANGELES,  CA 

DET  NORSKE  VERITAS  RESEARCH  AS  /  BERGAN,  VERITASVEIEN  1,  N-1322  HOVIK 

DTIC  /  ALEXANDRIA,  VA 

ENVIROPLEX  /  AUDIBERT,  HOUSTON,  TX 

FAU  /  BOCA  RATON,  FL 

FAU  /  REDDY,  BOCA  RATON,  FL 

GEORGE  WASHINGTON  UNIV  /  ENGRG  &  APP  SCI  SCHL  (FOX).  WASHINGTON,  DC 
GEOTECHNICAL  ENGRS,  INC  /  MURDOCK.  WINCHESTER.  MA 
GEOTEST  INSTRUMENT  CORP  /  BABENDREIER.  SILVER  SPRINGS.  MD 
GERWICK  INC  /  SAN  FRANCISCO.  CA 

GIANNOTTI  &  ASSOCIATES  OF  TEXAS  INC  /  GIANNOTTI,  VENTURA.  CA 

HEUZE  /  ALAMO,  CA 

HKS  INC  /  PAWTUCKET,  RI 

HQ  AFESC  /  TYNDALL  AFB,  FL 

IOWA  STATE  UNIV  /  AMES,  lA 

JOHN  HOPKINS  UNIV  /  COX,  BALTIMORE,  MD 

JOHN  HOPKINS  UNIV  /  LADE,  BALTIMORE,  MD 

KARAGOZIAN  AND  CASE  /  CRAWFORD,  GLENDALE,  CA 

LANTNAVFACENGCOM  /  CODE  411,  NORFOLK,  VA 


LAWRENCE  LIVERMORE  NATIONAL  LAB  /  MAKER,  LIVERMORE,  CA 

LAWRENCE  LIVERMORE  NATIONAL  LAB  /  MCCALLEN,  LIVERMORE,  CA 

LAWRENCE  TECH  UNIV  /  SOUTHFIELD,  ME 

LEHIGH  UNIV  /  BETHLEHAM,  PA 

LIN  OFFSHORE  ENGRG  /  SAN  FRANCISCO,  CA 

LOCKHEED  /  RSCH  LAB  (M.  JACOBY),  PALO  ALTO,  CA 

LOCKHEED  /  RSCH  LAB  (P  UNDERWOOD),  PALO  ALTO,  CA 

MAINE  MARITIME  ACADEMY  /  LIB,  CASTINE,  ME 

MARC  ANALYSIS  RSCH  CORP  /  HSU,  PALO  ALTO,  CA 

MCCLELLAND  ENGRS  /  HOUSTON,  TX 

MICHIGAN  UNIV  /  HOUGHTON,  MI 

MIT  /  R.V.  WHITMAN,  CAMBRIDGE,  MA 

MOBILE  R&D  CORP  /  DALLAS,  TX 

MONTANA  STATE  UNIV  /  PERKINS,  BOZEMAN,  MT 

NATL  ACADEMY  OF  SCIENCES  /  NRC,  DR.  CHUNG,  WASHINGTON,  DC 

NAVFACENGCOM  /  CODE  04 A3,  ALEXANDRIA,  VA 

NAVFACENGCOM  /  CODE  1002B,  ALEXANDRIA,  VA 

NAVFACENGCOM  /  CODE  163,  ALEXANDRIA,  VA 

NAVSHIP  /  CODE  245L,  FPO  AP, 

NCSC  /  PANAMA  CITY,  FL 

NEW  ENGLAND  MARINE  RESEARCH  LAB  /  LIB,  DUXBURY,  MA 
NFESC  ECDET  /  ESC61  (WU),  WASHINGTON,  DC 
NFESC  ECDET  /  ESC61  CECILIO,  WASHINGTON,  DC 
NORDA  /  BENNETT,  NSTL,  MS 

NORTH  CAROLINA  STATE  UNIV  /  RAHMAN,  RALEIGH,  NC 
NORTHROP  /  CHEN,  HAWTHORNE,  CA 
NORTHWESTERN  UNIV  /  LIU,  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  BAZANT,  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  CE  DEPT  (BELYTSCHKO),  EVANSTON,  IL 
NRL  /  VALENT,  STENNIS  SPACE  CENTER,  MS 
NSF  /  ASTILL,  ARLINGTON,  VA 

NSF  /  STRUC  &  BLDG  SYSTEMS  (KP  CHONG),  WASHINGTON,  DC 
NTH  /  LANGSETH,  N-7034  TRONDHEIM, 

NTH  /  MALO,  N-7034  TRONDHEIM, 

NUSC  /  LIB,  NEW  LONDON,  CT 

NY  MARITIME  COLLEGE  /  BRONX,  NY 

OCNR  /  CODE  10P4  (KOSTOFF),  ARLINGTON,  VA 

OCNR  /  CODE  1121  SILVA,  ARLINGTON,  VA 

ONR  /  CODE  1132SM,  ARLINGTON,  VA 

ONR  /  RAMBERG,  ARUNGTON,  VA 

OREGON  STATE  UNIV  /  CORVALLIS,  OR 

OREGON  STATE  UNIV  /  CORVALLIS,  OR 

OREGON  STATE  UNIV  /  CE  DEPT  (YIM),  CORVALLIS.  OR 

OREGON  STATE  UNIV  /  DEPT  OF  MECH  ENGRG  (SMITH),  CORVALLIS,  OR 

PENN  STATE  UNIV  /  LAB,  STATE  COLLEGE,  PA 

PMB  ENGRG  /  LUNDBERG,  SAN  FRANCISCO,  CA 

PORTLAND  STATE  UNIV  /  MIGLIORI,  PORTLAND,  OR 

PURDUE  UNIVERSITY  /  WEST  LAFAYETTE,  IN 

SAN  DIEGO  STATE  UNIV  /  SAN  DIEGO,  CA 

SCHWER  ENGR  &  CONSULTING  SERVICES  /  SCHWER,  BURR  RIDGE.  IL 
SCOPUS  TECHNOLOGY  INC  /  (B  NOUR-OMID),  EMERYVILLE,  CA 
SCOPUS  TECHNOLOGY  INC  /  (S.  NOUR-OMID),  EMERYVILLE,  CA 
SEAL  TEAM  6  /  NORFOLK,  VA 


SEATTLE  UNIV  /  SEATTLE,  WA 
SHANNON  AND  WILSON  INC  /  SEATTLE,  WA 
SJSU  /  VUKAZICH,  SAN  JOSE,  CA 

SOUTH  DAKOTA  SCHOOL  OF  MINES  AND  TECH  /  BANG,  RAPID  CITY,  SD 
SRI  INTL  /  ENGRG  MECH  DEPT  (GRANT),  MENLO  PARK,  CA 
SRI  INTL  /  ENGRG  MECH  DEPT  (SIMONS),  MENLO  PARK,  CA 
STANFORD  UNIV  /  APP  MECH  DIV  (HUGHES),  STANFORD,  CA 
STANFORD  UNIV  /  CE  DEPT  (PENSKY),  STANFORD,  CA 
STANFORD  UNIV  /  LAW,  STANFORD,  CA 
STATE  OF  CALIF  /  SACRAMENTO,  CA 

STRUCTURAL  ANALYSIS  PROGRAMS  INC  /  WILSON,  EL  CERRITO,  CA 

TEXAS  A&M  UNIV  /  COLLEGE  STATION,  TX 

TEXAS  A&M  UNIV  /  HERBICH,  COLLEGE  STATION,  TX 

TEXAS  A&M  UNIV  /  ROSCHKE,  COLLEGE  STATION,  TX 

THE  EARTH  TECH  CORP  /  ARULMOU,  IRVINE,  CA 

TU  DELFT  /  DE  BORST,  2600  GA  DELFT, 

TU  DELFT  /  VAN  MIER,  2600  GA  DELFT, 

TUFTS  UNIV  /  SANAYEI,  MEDFORD,  MA 
UCLA  /  JU,  LOS  ANGELES,  CA 

UCSD  /  SCRIPPS  INST  OF  OCEANOGRAPHY,  LA  JOLLA,  CA 
UNIV  OF  CAL  BERKELEY  /  FILIPPOU,  BERKELEY,  CA 
UNIV  OF  CAL  BERKELEY  /  GOVINDJEE,  BERKELEY,  CA 
UNTV  OF  CALIF  BERKELEY  /  BERKELEY,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (BAYO),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (BRUCH),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (LECKIE),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (MCMEEKING),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (TULIN),  SANTA  BARBARA,  CA 

UNTV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (HERRMANN),  DAVIS,  CA 

UNIV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (KUTTER),  DAVIS,  CA 

UNIV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (RAMEY),  DAVIS,  CA 

UNIV  OF  CALIFORNIA  DAVIS  /  CTR  FOR  GEOTECH  MODEL  GDRISS),  DAVIS,  CA 

UNTV  OF  COLORADO  /  CE  DEPT  (HON-YIM  KO),  BOULDER,  CO 

UNIV  OF  COLORADO  /  MECH  ENGRG  DEPT  (PARK),  BOULDER,  CO 

UNIV  OF  CONN  /  LEONARD,  STORRS,  CT 

UNIV  OF  CONN  /  LIBRARY,  GROTON,  CT 

UNIV  OF  DELAWARE  /  NEWARK,  DC 

UNIV  OF  HAWAII  /  HONOLULU.  HI 

UNTV  OF  HAWAII  /  KANEOHE  BAY,  HI 

UNIV  OF  HAWAII  /  HONOLULU,  HI 

UNTV  OF  HAW  An  /  HONOLULU,  HI 

UNIV  OF  ILLINOIS  /  CE  LAB  (PECKNOLD),  URBANA.  IL 

UNTV  OF  MICH  /  ANN  ARBOR.  ME 

UNIV  OF  N  CAROUNA  /  CE  DEPT  (GUPTA),  RALEIGH.  NC 
UNTV  OF  N  CAROLINA  /  CE  DEPT  (TUNG),  RALEIGH.  NC 
UNIV  OF  NY  /  BUFFALO,  NY 

UNIV  OF  RHODE  ISLAND  /  KOVACS,  KINGSTON,  RI 

UNIV  OF  RHODE  ISLAND  /  VEYERA,  KINGSTON,  RI 

UNIV  OF  TENNESSEE  /  KNOXVILLE,  TN 

UNIV  OF  TEXAS  /  JIRSA,  AUSTIN,  TX 

UNIV  OF  TEXAS  /  TASSOULAS,  AUSTIN,  TX 

UNIV  OF  WASHINGTON  /  MATTOCK,  SEATTLE,  WA 

UNIV  OF  WISCONSIN  /  M^WAUKEE,  WI 

UNTV  OF  WYOMING  /  CIVIL  ENGRG  DEPT,  LARAMIE,  WY 


USA  BELVOIR  /  FORT  BELVOIR,  VA 
USACOE  /  BRADLEY,  VICKSBURG,  MS 

VA  POLY  INST  AND  STATE  UNIV  /  MITCHELL,  BLACKSBURG,  VA 
VIRGINIA  INST  OF  MARINE  SCI  /  GLOUCESTER  POINT,  VA 
WEIDLINGER  ASSOCIATES  /  LEVINE,  LOS  ALTOS,  CA 
WEST  VIRGINIA  UNIV  /  BARBERO,  MORGANTOWN,  WV 
WEST  VIRGINIA  UNIV  /  KIGER,  MORGANTOWN,  WV 
WOODS  HOLE  OCEANOGRAPHIC  /  LIB,  WOODS  HOLE,  MA 
WOODWARD  CLYDE  CONSULTANTS  /  OAKLAND,  CA 
WORCHESTER  POLYTECH  /  SULLIVAN,  WORCESTER,  MA 


